{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Forward-Backward Filter]:  https://github.com/berndporr/iirj/issues/14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timeit\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import signal\n",
    "\n",
    "def sos_lowpass_filter(x, fs, high_cutoff_freq, order=12):\n",
    "    sos = signal.butter(order, high_cutoff_freq / (fs / 2.0), 'lowpass', output='sos')\n",
    "    y = signal.sosfiltfilt(sos, x)\n",
    "    return y\n",
    "\n",
    "def fft_ifft(x):\n",
    "    fft_size = len(x)\n",
    "    tmp = np.fft.fft(x, n=fft_size)\n",
    "    y = np.fft.ifft(tmp)\n",
    "    \n",
    "    # Remove imaginary part due to imprecisions\n",
    "    return np.real(y)\n",
    "\n",
    "def fft_lowpass_filter(x, fs, high_cutoff_freq):\n",
    "    fft_size = len(x)\n",
    "    \n",
    "    high_cutoff_idx = int(high_cutoff_freq * float(fft_size) / float(fs))\n",
    "    \n",
    "    tmp = np.fft.fft(x, n=fft_size)\n",
    "    tmp[high_cutoff_idx:fft_size-high_cutoff_idx+1] = 0\n",
    "    y = np.fft.ifft(tmp)\n",
    "    \n",
    "    # Remove imaginary part due to imprecisions\n",
    "    return np.real(y)\n",
    "\n",
    "\n",
    "# Timing comparison between sosfiltfilt and spectrum-based filtering\n",
    "    \n",
    "t = np.linspace(0.0, 60.0, num=60*44100, endpoint=False)\n",
    "x = np.random.randn(len(t))\n",
    "\n",
    "def sosfiltfilt_test():\n",
    "    return sos_lowpass_filter(x, 44100.0, 7500.0)\n",
    "\n",
    "def fft_filter_test():\n",
    "    return fft_lowpass_filter(x, 44100.0, 7500.0)\n",
    "\n",
    "t1 = timeit.timeit(\"sosfiltfilt_test()\", number=50, setup=\"from __main__ import sosfiltfilt_test\")\n",
    "print (\"IIR forward-backward time: %s\" % (t1,))\n",
    "t2 = timeit.timeit(\"fft_filter_test()\", number=50, setup=\"from __main__ import fft_filter_test\")\n",
    "print (\"Spectrum filtering time: %s\" % (t2,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample data output (lfilter doc sample data)\n",
    "\n",
    "t = np.linspace(-1, 1, 201)\n",
    "x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) +\n",
    "     0.1*np.sin(2*np.pi*1.25*t + 1) +\n",
    "     0.18*np.cos(2*np.pi*3.85*t))\n",
    "x = x + np.random.randn(len(t)) * 0.08\n",
    "\n",
    "x2 = fft_ifft(x)\n",
    "y1 = sos_lowpass_filter(x, 2.0, 0.05)\n",
    "y2 = fft_lowpass_filter(x, 2.0, 0.05)\n",
    "\n",
    "plt.figure(dpi=800)\n",
    "plt.plot(t, x, 'b', alpha=0.75, label='Signal')\n",
    "plt.plot(t, x2, 'k.', label='FFT -> IFFT')\n",
    "plt.plot(t, y1, 'm', label='Butterworth sosfiltfilt')\n",
    "plt.plot(t, y2, 'g', label='FFT -> Filter -> IFFT')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
